#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#define PI 3.14159265358979323846

/**
 * Sinc 插值函数（归一化 sinc: sin(pi*x)/(pi*x)）
 * @param t: 目标插值时间点（数组）
 * @param t_len: t 数组长度
 * @param x_n: 采样幅值数组
 * @param n_len: 采样点数量
 * @param T: 采样周期
 * @param result: 输出插值结果（长度 t_len）
 * @return: 0 成功，-1 失败
 */
int sinc_interp(const double *t, int t_len, const double *x_n, int n_len, double T, double *result) {
    if (t == NULL || x_n == NULL || result == NULL || t_len <= 0 || n_len <= 0 || T <= 0) {
        fprintf(stderr, "Invalid input parameters\n");
        return -1;
    }

    // 遍历每个目标插值点
    for (int i = 0; i < t_len; i++) {
        double sum = 0.0;
        // 遍历所有采样点，计算加权和
        for (int k = 0; k < n_len; k++) {
            double tau = (t[i] - k * T) / T;  // 归一化偏移量
            double sinc_val;
            if (fabs(tau) < 1e-12) {  // 避免除以零（tau=0 时 sinc(0)=1）
                sinc_val = 1.0;
            } else {
                sinc_val = sin(PI * tau) / (PI * tau);
            }
            sum += x_n[k] * sinc_val;
        }
        result[i] = sum;
    }
    return 0;
}

/**
 * 生成等间距时间轴（模拟 np.linspace）
 * @param start: 起始值
 * @param end: 结束值
 * @param num: 点数
 * @param out: 输出数组（需提前分配内存）
 */
void linspace(double start, double end, int num, double *out) {
    if (num <= 1) {
        out[0] = start;
        return;
    }
    double step = (end - start) / (num - 1);
    for (int i = 0; i < num; i++) {
        out[i] = start + i * step;
    }
}

/**
 * 调用 gnuplot 绘制图形
 * @param t_fine: 插值时间轴
 * @param x_fine: 插值结果
 * @param fine_len: 插值点数
 * @param x_n: 采样幅值
 * @param n_len: 采样点数
 * @param T: 采样周期
 */
void plot_result(const double *t_fine, const double *x_fine, int fine_len,
                 const double *x_n, int n_len, double T) {
    // 1. 生成临时数据文件
    FILE *data_file = fopen("sinc_interp_data.txt", "w");
    if (!data_file) {
        fprintf(stderr, "Failed to create data file\n");
        return;
    }

    // 写入插值曲线数据
    for (int i = 0; i < fine_len; i++) {
        fprintf(data_file, "%.6lf %.6lf\n", t_fine[i], x_fine[i]);
    }
    fprintf(data_file, "\n\n");  // 分隔符（gnuplot 多数据集）

    // 写入采样点数据
    for (int k = 0; k < n_len; k++) {
        fprintf(data_file, "%.6lf %.6lf\n", k * T, x_n[k]);
    }
    fclose(data_file);

    // 2. 构造 gnuplot 命令
    FILE *gnuplot = popen("gnuplot -persist", "w");
    if (!gnuplot) {
        fprintf(stderr, "Failed to start gnuplot (please install gnuplot first)\n");
        remove("sinc_interp_data.txt");  // 清理临时文件
        return;
    }

    // 设置中文字体（根据系统调整，Windows 可用 "SimHei"，Linux/Mac 可用 "WenQuanYi Micro Hei"）
    fprintf(gnuplot, "set encoding utf8\n");
    fprintf(gnuplot, "set terminal qt font 'WenQuanYi Micro Hei, 10'\n");
    fprintf(gnuplot, "set title 'Sinc 插值重建（周期 T=%.1lf）'\n", T);
    fprintf(gnuplot, "set xlabel '时间 t'\n");
    fprintf(gnuplot, "set ylabel '幅值'\n");
    fprintf(gnuplot, "set grid xtics ytics lw 1 lc rgb '#cccccc' dashtype '--'\n");
    fprintf(gnuplot, "set key left top\n");
    // 绘制插值曲线和采样点
    fprintf(gnuplot, "plot 'sinc_interp_data.txt' index 0 with lines lw 2 lc rgb 'steelblue' title 'Sinc 插值', \
                      '' index 1 with points pt 7 ps 1.5 lc rgb 'steelblue' pt 7 title '采样点'\n");

    pclose(gnuplot);
    remove("sinc_interp_data.txt");  // 清理临时文件
}

int main() {
    // ------------------- 1. 参数配置（对齐 Python 代码） -------------------
    const double T = 0.5;                  // 采样周期
    const double x_n[] = {7, 6, 7, -1, 7};// 采样幅值
    const int n_len = sizeof(x_n) / sizeof(x_n[0]);  // 采样点数量（5）
    const int fine_points = 2000;          // 细密采样点数
    const double t_start = 0.0;
    const double t_end = 12.0;             // 时间轴范围 0~12

    // ------------------- 2. 生成细密时间轴 -------------------
    double *t_fine = (double *)malloc(fine_points * sizeof(double));
    if (!t_fine) {
        fprintf(stderr, "Memory allocation failed for t_fine\n");
        return -1;
    }
    linspace(t_start, t_end, fine_points, t_fine);

    // ------------------- 3. 计算 Sinc 插值 -------------------
    double *x_fine = (double *)malloc(fine_points * sizeof(double));
    if (!x_fine) {
        fprintf(stderr, "Memory allocation failed for x_fine\n");
        free(t_fine);
        return -1;
    }
    int ret = sinc_interp(t_fine, fine_points, x_n, n_len, T, x_fine);
    if (ret != 0) {
        fprintf(stderr, "Sinc interpolation failed\n");
        free(t_fine);
        free(x_fine);
        return -1;
    }

    // ------------------- 4. 绘制图形 -------------------
    plot_result(t_fine, x_fine, fine_points, x_n, n_len, T);

    // ------------------- 5. 释放内存 -------------------
    free(t_fine);
    free(x_fine);

    printf("Sinc interpolation completed successfully\n");
    return 0;
}